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Time integration methods can be separated into two groups: explicit and 

implicit. Roughly speaking, methods which do not involve the solution of any 
algebraic equations are called explicit, whereas those that require the solution 
of equations are called implicit. 

The relative advantages and disadvantages of explicit and implicit methods 
are summarized in Fig. 1. It is interesting to observe that the positive attrib- 
utes of these two methods form complementary sets, so that if the positive 
attributes of the two methods can be combined into a single method, a truly 
powerful method would be achieved. 

An important point which is brought out in Fig. 1 is that whereas implicit 
methods are unconditionally stable for linear problems, stability does not imply 
accuracy and in fact the stability of implicit methods has often misled 
structural analysts into using time steps which yield very poor accuracy. 
Furthermore, no current time integration will undoubtedly be an important topic 
for future research. 


Relative merits of explicit and implicit integration methods 


Expl i ci t 

+ very simple and trouble free algorithm, complex phenomena easily included 
+ accuracy is assured if At stable for large systems 
+ no stiffness matrix necessary - saves storage 
conditionally stable, small At 

Implicit 

+ unconditionally stable, large At 

complex algorithm with low reliability in nonlinear situations 
accuracy can deteriorate 

Newton form has large core storage requirements 

Figure 1 
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The major trend of the past decade of research on time integration proce- 
dures has been hybridization methods so as to take advantage of the complementary 
nature of the positive attributes of explicit and implicit integration. The 
types of hybridization are indicated in Fig. 2. References for these methods are 
as follows: partitioning [1-7], operating splitting methods [8-11], semi- 
implicit methods [9-12], It should be noted that the distinction between semi- 
implicit methods and operator splitting methods is rather fuzzy; both groups of 
methods try to achieve unconditional stability through some modification of the 
evolution operator which either completely obviates the need for solving any 
equations or reduces the size of the system to be solved. 


Objective of current research in time integration: 

to exploit the advantages of implicit and explicit methods through 
hybridization (advantages of the two methods are complementary!) 

di rections: 

partitioning : different operators on different parts of the mesh 

semi-imp! icit methods : unconditionally stable methods that require no 
solution of equations or smaller systems 

operator splitting methods: split A to simplify solution - similar to 

semi -implicit 


Figure 2 
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The major shortcoming of operator splitting methods has been the rapid 
deterioration of their accuracy with increasing time step. For example, if we 
consider the Trujillo semi-implicit method, which is illustrated in Fig. 3, we 
find that as the Courant number increases the accuracy diminishes dramatically. 
In Reference [10] it is shown that the phase velocity in a one-dimensional mesh 
in the Trujillo method is such that the shorter waves essentially only advance 
one mesh length during a time step; thus, the effect of the semi-implicit inte- 
grator, as shown in Fig. 4, is to retard wave velocity so severely that regard- 
less of the size of the time step a quasi-Courant condition applies in that the 
numerical waves only traverse a single element in a time step. This distortional 
characteristic of semi-implicit methods has also been noted in the rigid-body 
modes by Park and Housner [12]. In Reference [12] several techniques for 
improving the accuracy of semi-implicit methods were developed, but we have not 
had time to check their effects independently. 


TRUJILLO SEMI-IMPLICIT 
(ref. 9) 

(M + it K) u n+1 = M u n + it s n+1 


let K = K l + Ky 


u 

(M + it K l ) u n+1 * = (M - it K u ) u n + it s n + 1 

(M + it K u ) u n+1 = (M - it K u ) u n+ ^ + it s n+1 

similar to 2 passes of Gauss - Seidel 

unconditionally stable 

accuracy? 

Figure 3 
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Avoidance of equation solution and unconditional stability can be achieved 
by rational Runge-Kutta methods [13], see Fig. 5. Again, the accuracy of these 
methods deteriorates rather quickly when the time step is much larger than the 
stability limit for the explicit methods. These methods seem to be most suited 
to parabolic systems. For structural mechanics , which involves a combination of 
hyperbolic and parabolic behavior, their lack of accuracy is generally unaccept- 
able. 


Rational Runqe Kutta 

« Vj I n+ 1 

« *2 + K(0 n + c 2 at vj) = f n+1 

~ = b l ~1 + b 2 ~2 
~n+l = ~n + at 

where e = 2(vj b) - (vj v ^ ) b 

unconditionally stable and second order accurate 

if c b x = 2, b 2 = -1 Hairer 1980 (ref. 13) 

no solution of equations if M diagonal 

x < 0 i f At is too large 

partitioned Rational Runge Kutta methods, Liu et al . 

I JNME , 1581-1597, (1984), 1984 (ref. 14) 


Figure 5 
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A very novel operating splitting method, which exploits the unique features 
of the finite-element assembly operation, has recently been developed by Hughes 
and coworkers [15] . This method only required conversion of the element matri- 
ces, so while the method does not completely avoid the solution of equations as 
in semi-implicit methods, the size of equations to be solved is reduced substan- 
tially, see Fig. 6. Hughes and coworkers make a very compelling argument that 
this type of method will prove particularly beneficial in three-dimensional 
applications . 

We have tested an early version of the method in both parabolic systems and 
elastic-plastic structural mechanics problems. In comparing the method with a 
conjugate gradient method, we found that the element-by-element and conjugate 
gradient methods were of comparable speed. When used with large time steps in 
structural dynamics problems, we were not able to achieve reasonable accuracy 
unless we made a large number of sweeps during each time step. On the other 
hand, we found the method to be very useful in crash-type problems in conjunction 
with explicit techniques. As a deforming structure becomes mostly plastic, it 
becomes possible to increase the explicit time step quite a bit if the eleraent- 
by-element procedure is used to stabilize elements which unload into the elastic 
regime. This would detract somewhat from the accuracy if it occurs in many ele- 
ment. However, in general, phase accuracy is not an overriding concern in crash- 
type problems, so that the potential of these methods for stabilizing explicit 
methods is worth investigating. We have not yet tried the later versions of the 
eleraent-by-element technique which are reported to be substantially more accu- 
rate. Reference [16] reports a procedure which reduces the computational effort 
required in solving the element equations by as much as an order of magnitude. 

ELEMENT-BY-ELEMENT OPERATOR SPLIT 


Hughes, Levlt, Wlnget A5CE J. Eng. Mech. Oiv. April 1983 (ref. 8) 
Comp. Meth. Appl. Mech. Eng. 1983 (ref. 15) 
Ortiz, Pinsky and Taylor (ref. 17) 

Recall implicit Eqns. 

(M + At JC) u n+l » 7 (1) 

(£ ♦ At K) u * f a M' 1 f (superscript dropped) 

( 2 ) 

Approximation 


(I ♦ At M' 1 l JC ) * J| 
e e 


G 

-e 


So (2) becomes 


n 5 e id 5 (5l £2 •••• £ne^ £ * l 

Procedure 


-CO] * - f S[.] ' U Ce-13 

ONLY ELEMENT MATRICES NEEO TO BE INVERTED! 

Figure 6 
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For problems with different time scales such as the space-structure 
deployment problem, where high-frequency impacts occur in conjunction with low- 
frequency rigid-body motions, the partitioned methods are quite promising. 
Partitioned methods are defined as those which employ different time steps or 
different integrators in different parts of the mesh. During an input, it would 
be desirable to use different time steps in the vicinity of the impact in solving 
a large-scale structure problem. By doing this, accuracy could be retained in 
all parts of the mesh without engendering large expense. The potential of these 
methods is indicated in Fig. 7. 


At a = 0.14 MSEC At ST = 0.04 MSEC 

RUNNING TIME FOR 60 MSEC SIMULATION 



FOR 60 SEC SEISMIC SIMULATION 
RUNNING TIME E - I: 4.2 X 10 5 SEC 


Figure 7 
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Considerable progress has recently been achieved in mesh partitions with 
different time steps, see References 5 to 12; 14 to 16, and 18 to 20. Basically, 
two types of mixed time partitions have been involved: element partitions and 
nodal partitions. The algorithm for nodal partition is shown in Fig. 8. Nodal 
partitions appear to provide the best accuracy, but their analysis has been 
impeded by the fact that the amplification matrix is not symmetric. 



• amplification matrix is not symmetric 

Figure 8 
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An elemental partition is shown in Fig. 9. Element partitions are 
associated with symmetric amplification matrices and in Ref. [20] a proof of 
sufficient conditions or stability has been given for a first-order, linear 
system with different time steps. The proof applies to both explicit and 
implicit integrators. 


Subcyclinq with Elemental Partition 


u 


n +1 


at (if 1 )!/ 



1 st order system i.e. heat conduction 
di ffusion 
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O' 
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3 4 


elements 1 - with At 

elements 2 & 3 - with 2At 


computations in cycle 


update 

update 

update 

update 

update 

update 


u l> u 2 

u, , u 0 i i sometimes deleted 

1 L£j 

f 0 i = 1 to 3 

Up u 2> u 3> u 4 

f © 


amplification matrix eqns are symmetric 


Figure 9 
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Partitioned implicit methods are particularly well-suited to iterative 
solvers. Whereas for Newton-type solvers, several different triangulations have 
to be stored for mixed time integration, this is not necessary for iterative 
solvers. To illustrate the nature of the solutions which can be obtained from 
these methods the results of the thermal transient problem in Fig. 10 are shown 
in Figs. 11 and 12, An interesting observation from Fig. 12 is that when the 
time step ratio is extremely large (32:1 in case 2), stability is maintained but 
large errors develop. It has become clear that methods of this type must use a 
smooth transition of time steps from the smallest time step to the largest time 
step. Thus, an important ingredient in any mixed time integration procedure is a 
strategy which automatically selects the time steps within the different regimes 
according to accuracy requirements and provides smooth transitions of time steps 
between regions where very large time steps can be used and those where very 
small time steps can be used. 




Figure 10 






The potential of these methods even in two-dimensional problems is quite 
tremendous, as evidenced by the comparisons shown in Fig, 13. This is a two- 
dimensional heat conduction problem with a large range of thermal conductivities. 
As can be seen from the comparison, savings of a factor of 2 to 5 can be achieved 
even in moderately sized two-dimensional problems. These types of savings have 
important implications in a computer-aided engineering environment, where the 
analysis of a new concept must be achieved in a reasonable amount of time if the 
design process is to be interactive. 

These mixed-time integration procedures are in many ways still in their 
infancy. The applications to nonlinear problems and contact-impact problems will 
probably require special strategies in order to exploit these methods to their 
fullest advantage. It would also be desirable to develop stability analyses in 
the linear regime for second-order systems, such as the equations of motion, and 
for nodal partitions. 

This class of methods, when combined with iterative solvers, would be 
uniquely suited for parallel architecture computers. In principle, each sub- 
domain with a particular time step could be treated by a different CPU. Data 
transfer between subdomains would only be necessary for interface data. 


Storage and Running Time Comparisons 


Storage 

Method I-4E-8E-8I I I-E-2E-2I I 

At - 1 At - 1 At = 4 At = 4 


nonzero 
terms in K 

3288 

33771 

6089 

33771 

average 
semi bandwidth 

6 

69 

12 

69 

solution 

time 

CPU-s 

90. 

516. 

70. 

139. 


Note: • problem is linear; 

e 8 x 50 mesh is numbered for large bandwidth to simulate 
3D problems. 


Figure 13 
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